{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Introduction to ITK Segmentation in SimpleITK Notebooks\n",
    "\n",
    "<b>Goal</b>: To become familiar with basic segmentation algorithms available in ITK, and interactively explore their parameter space.\n",
    "\n",
    "Image segmentation filters process an image to partition it into (hopefully) meaningful regions. The output is commonly an image of integers where each integer can represent an object. The value 0 is commonly used for the background, and 1 ( sometimes 255) for a foreground object.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from ipywidgets import interact, FloatSlider\n",
    "\n",
    "import SimpleITK as sitk\n",
    "\n",
    "# Download data to work on\n",
    "%run update_path_to_download_script\n",
    "from downloaddata import fetch_data as fdata\n",
    "from myshow import myshow, myshow3d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "img_T1 = sitk.ReadImage(fdata(\"nac-hncma-atlas2013-Slicer4Version/Data/A1_grayT1.nrrd\"))\n",
    "img_T2 = sitk.ReadImage(fdata(\"nac-hncma-atlas2013-Slicer4Version/Data/A1_grayT2.nrrd\"))\n",
    "\n",
    "# To visualize the labels image in RGB with needs a image with 0-255 range\n",
    "img_T1_255 = sitk.Cast(sitk.RescaleIntensity(img_T1), sitk.sitkUInt8)\n",
    "img_T2_255 = sitk.Cast(sitk.RescaleIntensity(img_T2), sitk.sitkUInt8)\n",
    "\n",
    "myshow3d(img_T1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Thresholding\n",
    "\n",
    "Thresholding is the most basic form of segmentation. It simply labels the pixels of an image based on the intensity range without respect to geometry or connectivity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "seg = img_T1>200\n",
    "myshow(sitk.LabelOverlay(img_T1_255, seg), \"Basic Thresholding\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "seg = sitk.BinaryThreshold(img_T1, lowerThreshold=100, upperThreshold=400, insideValue=1, outsideValue=0)\n",
    "myshow(sitk.LabelOverlay(img_T1_255, seg), \"Binary Thresholding\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ITK has a number of histogram based automatic thresholding filters including Huang, MaximumEntropy, Triangle, and the popular Otsu's method. These methods create a histogram then use a heuristic to determine a threshold value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "otsu_filter = sitk.OtsuThresholdImageFilter()\n",
    "otsu_filter.SetInsideValue(0)\n",
    "otsu_filter.SetOutsideValue(1)\n",
    "seg = otsu_filter.Execute(img_T1)\n",
    "myshow(sitk.LabelOverlay(img_T1_255, seg), \"Otsu Thresholding\")\n",
    "\n",
    "print(otsu_filter.GetThreshold() )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Region Growing Segmentation\n",
    "\n",
    "The first step of improvement upon the naive thresholding is a class of algorithms called region growing. This includes:\n",
    "<ul>\n",
    "  <li><a href=\"http://www.itk.org/Doxygen/html/classitk_1_1ConnectedThresholdImageFilter.html\">ConnectedThreshold</a></li>\n",
    "  <li><a href=\"http://www.itk.org/Doxygen/html/classitk_1_1ConfidenceConnectedImageFilter.html\">ConfidenceConnected</a></li>\n",
    "  <li><a href=\"http://www.itk.org/Doxygen/html/classitk_1_1VectorConfidenceConnectedImageFilter.html\">VectorConfidenceConnected</a></li>\n",
    "  <li><a href=\"http://www.itk.org/Doxygen/html/classitk_1_1NeighborhoodConnectedImageFilter.html\">NeighborhoodConnected</a></li>\n",
    "</ul>\n",
    "\n",
    "Earlier we used 3D Slicer to determine that index: (132,142,96) was a good seed for the left lateral ventricle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "seed = (132,142,96)\n",
    "seg = sitk.Image(img_T1.GetSize(), sitk.sitkUInt8)\n",
    "seg.CopyInformation(img_T1)\n",
    "seg[seed] = 1\n",
    "seg = sitk.BinaryDilate(seg, 3)\n",
    "myshow(sitk.LabelOverlay(img_T1_255, seg), \"Initial Seed\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "seg = sitk.ConnectedThreshold(img_T1, seedList=[seed], lower=100, upper=190)\n",
    "\n",
    "myshow(sitk.LabelOverlay(img_T1_255, seg), \"Connected Threshold\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Improving upon this is the ConfidenceConnected filter, which uses the initial seed or current segmentation to estimate the threshold range."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "seg = sitk.ConfidenceConnected(img_T1, seedList=[seed],\n",
    "                                   numberOfIterations=1,\n",
    "                                   multiplier=2.5,\n",
    "                                   initialNeighborhoodRadius=1,\n",
    "                                   replaceValue=1)\n",
    "\n",
    "myshow(sitk.LabelOverlay(img_T1_255, seg), \"ConfidenceConnected\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "img_multi = sitk.Compose(img_T1, img_T2)\n",
    "seg = sitk.VectorConfidenceConnected(img_multi, seedList=[seed],\n",
    "                                             numberOfIterations=1,\n",
    "                                             multiplier=2.5,\n",
    "                                             initialNeighborhoodRadius=1)\n",
    "myshow(sitk.LabelOverlay(img_T2_255, seg))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fast Marching Segmentation\n",
    "\n",
    "The FastMarchingImageFilter implements a fast marching solution to a simple level set evolution problem (eikonal equation). In this example, the speed term used in the differential equation is provided in the form of an image. The speed image is based on the gradient magnitude and mapped with the bounded reciprocal $1/(1+x)$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "seed = (132,142,96)\n",
    "feature_img = sitk.GradientMagnitudeRecursiveGaussian(img_T1, sigma=.5)\n",
    "speed_img = sitk.BoundedReciprocal(feature_img) # This is parameter free unlike the Sigmoid\n",
    "myshow(speed_img)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output of the FastMarchingImageFilter is a <b>time-crossing map</b> that indicates, for each pixel, how much time it would take for the front to arrive at the pixel location."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fm_filter = sitk.FastMarchingBaseImageFilter()\n",
    "fm_filter.SetTrialPoints([seed])\n",
    "fm_filter.SetStoppingValue(1000)\n",
    "fm_img = fm_filter.Execute(speed_img)\n",
    "myshow(sitk.Threshold(fm_img,\n",
    "                    lower=0.0,\n",
    "                    upper=fm_filter.GetStoppingValue(),\n",
    "                    outsideValue=fm_filter.GetStoppingValue()+1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def fm_callback(img, time, z):\n",
    "    seg = img<time;\n",
    "    myshow(sitk.LabelOverlay(img_T1_255[:,:,z], seg[:,:,z]))\n",
    "           \n",
    "interact( lambda **kwargs: fm_callback(fm_img, **kwargs),\n",
    "            time=FloatSlider(min=0.05, max=1000.0, step=0.05, value=100.0),\n",
    "            z=(0,fm_img.GetSize()[2]-1))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Level-Set Segmentation\n",
    "\n",
    "There are a variety of level-set based segmentation filter available in ITK:\n",
    "<ul>\n",
    "<li><a href=\"http://www.itk.org/Doxygen/html/classitk_1_1GeodesicActiveContourLevelSetImageFilter.html\">GeodesicActiveContour</a></li>\n",
    "<li><a href=\"http://www.itk.org/Doxygen/html/classitk_1_1ShapeDetectionLevelSetImageFilter.html\">ShapeDetection</a></li>\n",
    "<li><a href=\"http://www.itk.org/Doxygen/html/classitk_1_1ThresholdSegmentationLevelSetImageFilter.html\">ThresholdSegmentation</a></li>\n",
    "<li><a href=\"http://www.itk.org/Doxygen/html/classitk_1_1LaplacianSegmentationLevelSetImageFilter.html\">LaplacianSegmentation</a></li>\n",
    "<li><a href=\"http://www.itk.org/Doxygen/html/classitk_1_1ScalarChanAndVeseDenseLevelSetImageFilter.html\">ScalarChanAndVese</a></li>\n",
    "</ul>\n",
    "\n",
    "There is also a <a href=\"http://www.itk.org/Doxygen/html/group__ITKLevelSetsv4.html\">modular Level-set framework</a> which allows composition of terms and easy extension in C++.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we create a label image from our seed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "seed = (132,142,96)\n",
    "\n",
    "seg = sitk.Image(img_T1.GetSize(), sitk.sitkUInt8)\n",
    "seg.CopyInformation(img_T1)\n",
    "seg[seed] = 1\n",
    "seg = sitk.BinaryDilate(seg, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the seed to estimate a reasonable threshold range."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "stats = sitk.LabelStatisticsImageFilter()\n",
    "stats.Execute(img_T1, seg)\n",
    "\n",
    "factor = 3.5\n",
    "lower_threshold = stats.GetMean(1)-factor*stats.GetSigma(1)\n",
    "upper_threshold = stats.GetMean(1)+factor*stats.GetSigma(1)\n",
    "print(lower_threshold,upper_threshold)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "init_ls = sitk.SignedMaurerDistanceMap(seg, insideIsPositive=True, useImageSpacing=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lsFilter = sitk.ThresholdSegmentationLevelSetImageFilter()\n",
    "lsFilter.SetLowerThreshold(lower_threshold)\n",
    "lsFilter.SetUpperThreshold(upper_threshold)\n",
    "lsFilter.SetMaximumRMSError(0.02)\n",
    "lsFilter.SetNumberOfIterations(1000)\n",
    "lsFilter.SetCurvatureScaling(.5)\n",
    "lsFilter.SetPropagationScaling(1)\n",
    "lsFilter.ReverseExpansionDirectionOn()\n",
    "ls = lsFilter.Execute(init_ls, sitk.Cast(img_T1, sitk.sitkFloat32))\n",
    "print(lsFilter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "myshow(sitk.LabelOverlay(img_T1_255, ls>0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
